A discrete inhomogeneous model for the yeast cell cycle 
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We study the robustness and stability of the yeast cell regulatory network by using a general 
. . . inhomogeneous discrete model. We find that inhomogeneity, on average, enhances the stability of 

00 ' the biggest attractor of the dynamics and that the large size of the basin of attraction is robust 

\ against changes in the parameters of inhomogeneity. We find that the most frequent orbit, which 

. represents the cell-cycle pathway, has a better biological meaning than the one exhibited by the 

homogeneous model. 
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I. INTRODUCTION 



^ ' The eukaryotic cell exhibits a common process of division into two daugther cells. This process consists of four 
^ \ phases (i) Gl phase, in which the cell grows; (ii) S phase, in which the DNA is replicated; (iii) G2 phase, that 
^ ■ is a temporally gap between the S phase and the next one; (iv) M phase, in which the cell divides itself in two cells. 
Q I In Gl phase the cell cycle rests in a stationary state until the cell size reachs a critical value and the cell receive 
external signals which allow it to go on the cycle. The cell-cycle regulation machinery, which controls the growth and 
division processes, is known for the budding yeast in detail, Saccharomyces cerevisiae [3, B . In order to understand 
O^, the budding yeast regulation, several models have been proposed and discussed [1, H, U [a B S B B [13, El • Chen 
et al. [2] converted the regulation mechanism into a set of differential equations with empirical parameters. Their 
i kinetic model has taken in account several physiological, biochemistry and genetical details. Recently, Li et al. Q 
J> ■ introduced a simple Boolean dynamical model to investigate the stability and robustness features of the regulatory 
Tij" ' network. They found that the cell-cycle network is stable and robust. A stochastic version of this model, controlled by 
C ■ ■ a temperature like parameter, was latter studied p^]. The authors found that the stationary state and the cell-cycle 
' pathway are stable for a wide range of the temperature parameter values. Other aspects related to the checkpoint 
^; . efficiency of the Li et al. model were also considered p^. 

' In this work we consider the inhomogeneous version of the yeast cell-cycle introduced by Li et al. ^3,]. The cell 
cycle is represented by a regulatory network and the dynamics is modeled as a simple discrete dynamical system. The 
dynamics is constrained by the network topology by means of some parameters (coefficients) related to each network 
link. The link between nodes i and j determines the value of parameter Kij, which is present in the time evolution 
^ • rule. If the coefficient is positive, the link represents an activation and it will be noted as Kij = a{i,j) {a{i,j) > 0). 
k>( ' On the other hand, a negative coefficient Kij = —b{i,j) {b{i,j) > 0) represents an inactivation link. In the paper 
''jjj of Li et al. @, the authors studied basically the homogeneous model a{i,j) — b{i,j) = 1 for the nodes with a link 
between them. At most they have considered only two kinds of intensity a{i,j) = Ur and b(i,j) — ag. However, the 
interactions are important for the dynamics and are related to the constant rates of the kinetic equations, implying 
that different contributions to activation or inactivation may be present. We consider the most general inhomogeneous 
model of the yeast cell-cycle. Since that several different inhomogeneities represent the same dynamical model, we 
eliminate all possible degeneracy by constructing a minimal set of parameters (coefficients) . Such set represents all 
kind of inhomogeneity and is very large. We find that the big basin of attraction corresponding to the stationary 
state is still robust to changes in the coefficients and that in the minimal set of coefficients a new orbit corresponding 
to the cell cycle appears more frequently having a more feasible biological significance. This means that this orbit is 
more robust against change in coefficients than the orbit of the homogeneous model. Moreover, the mean basin size 
of the global attractor is bigger than for the homogeneous case, when this more frequent orbit is present. 



II. THE DYNAMICAL MODEL 



Our model is based on the deterministic Boolean model of Li et al. 's'l for the budding yeast cell-cycle regulation, 
which is represented by a regulatory network [1] Each one of the 11 nodes, which represents proteins or protein 
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FIG. 1: The yeast cell cycle regulatory network Sj. The nodes are identified by the proteins or complex proteins and by the 
number used in the definition of the dynamical system. The continuous arrows mean activation and the dashed ones mean 
inhibition. 



complexes, is represented by a variable S that takes the values (the protein state is inactive) or 1 (the protein state 
is active). The configuration of the system at time t is described by a vector S{t) = {Si{t), 5'2(i), . . . , Sii{t)) that 
represents the state of the following proteins or complex proteins: (Cln3, MBF, SBF, Clnl-2, Cdhl, Swi5, Cdc20, 
Clb5-6, Sicl, Clbl-2, Mcml). A configuration can be expressed in a short fashion by an integer I if we define that 
/ = J2i=i ■ -f'o^ example, if only Cdhl and Sicl are active {S5 — 1 and 5*9 = 1), the related configuration 
S = (0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0) corresponds to I{S) = 272. 

The state vector S describes the cell state in a specific time. If one wants to capture the time evolution of the cell 
states, one needs to address a dynamical model, in which there is a time evolution rule. The configuration at time 
t + 1, S{t + 1) is related to the previous one S{t) by the relationships 

r 1, ifEjK.,,s,{t)>o 
s.{t + i) = { 0, ifE,^.,.^.W <o (1) 
[ s^it), ifE,^.,j-^j(i) = , 

where Kij is the interaction between nodes i and j and i = 1, 2, . . . , 11. If there is no link between the nodes, we have 
that Kij = 0. Let us set our notation. When this link represents an activation, Kij = +a{i,j). Here a{i,j) is the 
intensity of such activation. When the interaction represents an inhibition we have that Ki j = —b{i,j), with b{i,j) 
being the inhibition intensity. Note that nodes 1, 4, 6, 7 and 11 have also a time-delayed self-degradation mechanism. 
When S'j(t) = 1 (i = 1, 4, 6, 7, 11) and J^j ^ij^j = from time t to t + td, we wih have that Si{t + td) = 0. From 
now on we will consider td ^ 1. This evolution rule can be set in a more compact form: 

r 1, ifa;>0 

S^{t + l)=F, ^ . K,,,S,{t)] , where F,{x) = <^ 0, if a; < (2) 
^ ^ [ 5,(0, if a: = . 

If we name Vl as the entire state space vector, which is the space that contains the 2^^ — 2048 possible vectors S{t), 
and if we name F the map defined by S{t + 1) = F{S{t)) — {Fi{Si{t)), . . . , Fii{Sii{t))), we can define a dynamical 
systems as the pair {il,F), em que F : ^ il. 

From Fig. [T] we can obtain the dynamical equations. For further use, let us write these equations as the trivial 
ones, namely 

Si{t + 1) - , 

Si{t + 1) = Fi{a{4,3)S3{t)) , 

S7{t + 1) = F7{a{7,10)Sw{t)+a{7,n)Siiit)) , (3) 
Sii{t + 1) - Fniain,8)Ssit)+a{ll,10)SMt)) , 
those involving one positive and one negative links 



S2{t + l) = F2{a{2,l)Si{t)-b{2,10)Sio{t)) 



(4) 
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F3{a{3,l)Si{t)-b{3,10)Sio{t)) 



(5) 



the ones involving one positive and two negative links, and the symmetrical case. 



Ssit + l) 



FcMQ, 7)S7{t) + a(6, ll)5ii(t) - b{6, lO)Sio{t)) 
Fsia{8,2)S2{t) - 6(8, 7)57(t) -6(8,9)59(0) , 



(6) 
(7) 



one equation involving three negative and one positive links 



S^it + 1) = F5(a(5, 7)Sj{t) - 6(5, 4)^4(0 - 6(5, 8)Ss{t) - 6(5, 10)5io(0) 



(8) 



and, finally, those involving three negative and two positive parameters 




of results. Now we will study the general inhomogeneous case. Set new values to the coefficients means change the 
dynamical system definition. But will these value changes actually define new dynamical systems? Let C/ be a subset 
of r2. If F([/) — F'{U) for all subset U, the dynamical systems (ft, F) and (fi, F') are the same. In the inhomogeneous 
model, the dynamical system is defined by the set C — {Kij;Q < i < 11,0 < j < 11}. Is that possible that a class 
[C], which can contain more than one set of coefhcients C, defines the same dynamical system? If this happens the 
non redundants inhomogeneities can be set by choosing only one element of each class [C] . 

In order to find this minimal set we must consider tow points: (i) the rule for time evolution takes into account only 
the sign of the sum J2j ^'^d (ii) the dynamical systems A and B are the same if, and only if, SA{t + 1) = 

Ssit + 1) for all possible common initial condition S{t) in the state space. 

We will illustrate how to find the minimal set explicitly for the dynamic of node 7, which is in equation [3l Note 
that both terms in the sum have the same sign. Whatever the values assumed by the pair C = {a(7, 10), a(7, 11)}, 
all possible initial state combinations (S'io(i), Sii{t)), namely D = {(0, 0), (1, 0), (0, 1), (1, 1)}, will result on the same 
case (positive, negative or zero) for the sum KijSj{t), namely S = {0, +, +, +}. But only the sign of the sum is 
important for the evaluation of 6*7 (t + 1), and we have that all the four possible values of the pair C give the same 
map between S and D. If two dynamical system have the same map between S and D they are identical. So we 
have all the four possible pairs defining the same dynamical systems what means that they take part of the same 
class [C]. We only need to take a representative pair of this class, being the simplest way to choose C = (1, 1). For 
all the other nodes the procedure is the same, but we use computer programs to make handle the calculations. The 
equation [3] have only one class, with all coefficients set equal one. The equations [4] and [5] have 3 classes, represented 
by {(1, 1), (2, 1), (1, 2)}. The equation [6] and [7] have 11 classes, represented by 



The equation [8] has 67 classes and the equations [9] and [TO] have 3265 classes. The minimal set that accounts for all 
the possible inhomogeneities has 3^ • 11^ ■ 67 ■ 3265^ = 7, 77 x 10^""^ elements. This huge quantity means that we must 
use a numerical simulation in order to sample this huge set. 



First let us remind briefly the main results obtained by Li et al. [3| for the homogeneous model with td — 1- From the 
2, 048 initial conditions, 1, 764 converge to the fixed point Ii = 272 representing the biological Gi stationary state. The 
other initial conditions converge to six fixed points {I2 = 0, I3 — 12, I4 — 16, I5 — 256, /g = 258 and Ij = 274). The 
pathway of the cell-cycle network is started by perturbing the Gi stationary state: Cln3 is turned on. Then the state of 
the model returns to the Gi fixed point after 12 time steps. The evolution of the proteins (or protein complexes), the so- 
called pathway of the cell-cycle network Pathhomog = [273, 278, 286, 14, 142, 1678, 1736, 1632, 1888, 1376, 368, 304, 272] , 
follows the biological cell-cycle sequence. 

We study the inhomogeneous model by randomly choosing the coefficients from the minimal set of interactions. 
In a typical numerical study we have 10^ samples. For each set of coefficients we study the evolution from each 
one of the 2,048 initial conditions and determine the fixed points (or the rare period-2 and period-3 cycles) with its 
attraction basin. The Gi fixed point 272 is always present. In our simulations we find that its minimal basin size 
is BSmin{'272) — 11. It occurs only in 0.34% of the minimal interaction set. More rare, the maximal basin size 



{(1,1,1), (1,1,2), (1,1,3), (1,2,1), (1,2,2), (1,3,2 



:), (2, 1,1), (2, 1,2), (2, 2,1), (2, 2, 3), (3, 1,2)}. 



III. RESULTS FOR THE INHOMOGENEOUS MODEL 
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Time Cln3 Mbf SBF Clnl-2 Cdhl Swi5 Cdc20 Clb5-6 Sicl Clbl-2 Mcml Phase I 
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TABLE I: Orbit with initial configuration being the perturbed stationary state. The column "Time" corresponds to the time 
steps; the columns identified by the nodes correspond to the time evolution of this nodes; the column "Phase" corresponds to 
the cell cycle phases; the column "I" correspond to the integer notation for the vector states. For Cdhl in times 8 and 9 the 
notation 0/1 means that the left state corresponds to Pathhomog and the right one to Pathfreq- 



BSmax{'272) — 1999 found in our simulations, appears in 0.23% of the interactions cases. In this case, generally we 
have only 5 fixed points, instead of 7 found in the homogeneous model. The mean basin size of the Gi fixed point, 
BSmean{272) = 1622.9, is smaller than that of the Li et al. model {BShomog{212) = 1764). However, 65.3% of the 
interaction's minimal set have BS{272) > 1764. 

We study also the pathway of the cell-cycle network. For each sample of coefficients, the temporal evolution of 
the pathway begins with the perturbed state 273, which is the stationary Gl state (272) with the Cln3 activation. 
Then we determine the pathway and its time step. In 10^ samples we find 336 different pathways. The pathway 
Pathfreq = [273, 278, 286, 14, 142, 1678, 1736, 1632, 1904, 1392, 368, 304, 272] with 12 time steps is the most frequent. 
It appears in 10.9% of the minimal interaction set, while the pathway of the Li et al. model appears only in 3.6% of the 
cases. If we evaluate the average basin size of the Gi fixed point in the ensemble of the cases having the most frequent 
pathway, we find that BSfreq{272) — 1866.7. Note that this average value is large than BShomog{272) = 1764. The 
only difference between Pathfreq and Pathhomog is that in the phase M the cyclin complex Cdhl is turned on at time 
step 9 in Pathfreq instead of time step 11 in Pathhomog, as it is shown in table [H Is this difference meaningful!? 

The exit of phase M is controlled by the complex APC, which is activated mainly by Cdc20 and Cdhl. The protein 
Cdc20 actives APC that begins the degradation of Clb2 during the transition metaphase-anaphase. The second 
phase of the Clb2 degradation occurs by the linking of Cdhl with APC, during the telophase. Therefore, the Clb2 
degradation needs a sequential action [Ij]: (i) Cdc20 activating APC and (ii) Cdhl linking to APC. Although the 
degradation by APC-Cdc20 is enough to the exit of phase M, the degradation by APC-Cdhl is essential to keep the 
Gl phase stable [l3|- In the homogeneous model Cdhl is turned on simultaneously with the turning off of Clb2, 
implying that the second degradation mechanism is absent. On the other hand, this second mechanism is present in 
Pathfreq bccausc Cdhl and Clb2 are both turned on for two time steps. We can conclude that the early activation 
of Cdhl is more coherent with the Clb2 degradation mechanism. 

This results show that not just the big basin of attraction of the stationary state is indeed robust against changes 
in the coefficients, but with inhomogeneities this basin of attraction gets bigger: 65, 5% of the interaction's minimal 
set have basin size bigger than 1764, with a maximum value of 1999. 

IV. CONCLUSIONS 

We found that the big basin of attraction is indeed robust against change in the parameter Ki^: 65, 3% of the min- 
imal interaction set have BS{272) > 1764 and the maximum size being 1999. Considering the orbit of the perturbed 
stationary state 273 we found that the most frequent orbit (Pfreq) is not the one exhibited by the homogeneous model: 
the orbit Pfreq appears in 10, 9% of the minimum coefficient set while Phomog appears only in 3, 6%, that is, this orbit 
is more robust against changes on the coefficients. Besides that, the subset of the dynamical systems exhibiting Pfreq 
has a mean basin size BSfreq = 1866,7 that is bigger than BShomog — 1764. Not just more robust, but Pfreq has a 
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better biological significance as we have shown by means of the earlier activation of Cdcl. 
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